{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we'e going to actually run a short simulation with OpenMM\n",
    "saving the results to disk with MDTraj's HDF5 reporter\n",
    "\n",
    "Running this example calculation on your machine requires\n",
    "having OpenMM installed. OpenMM can be downloaded and installed from\n",
    "https://simtk.org/home/openmm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets import some things we're going to need from mdtraj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import mdtraj\n",
    "import mdtraj.reporters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And a few things from OpenMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from simtk import unit\n",
    "import simtk.openmm as mm\n",
    "from simtk.openmm import app"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's find a PDB for alanine dipeptide, the system we'll\n",
    "be simulating. We've provided a PDB file under `data/` in this directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdb = mdtraj.load('data/native.pdb')\n",
    "topology = pdb.topology.to_openmm()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets use the `amber99sb-ildn` forcefield with implicit solvent\n",
    "and a Langevin integrator. This is a relatively \"standard\" OpenMM\n",
    "protocol for setting up a system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forcefield = app.ForceField('amber99sbildn.xml', 'amber99_obc.xml')\n",
    "system = forcefield.createSystem(topology, nonbondedMethod=app.CutoffNonPeriodic)\n",
    "integrator = mm.LangevinIntegrator(330*unit.kelvin, 1.0/unit.picoseconds, 2.0*unit.femtoseconds)\n",
    "simulation = app.Simulation(topology, system, integrator)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the initial positions to the \"first frame\" of the PDB\n",
    "file (it only has one frame). Note that while the pdb is an mdtraj\n",
    "trajectory passing in its positions to OpenMM is just fine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation.context.setPositions(pdb.xyz[0])\n",
    "simulation.context.setVelocitiesToTemperature(330*unit.kelvin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use one of the OpenMM reporters that mdtraj provides. This is\n",
    "the hdf5 reporter, which saves all kinds of information, including\n",
    "the topology, positions, energies, etc to disk. To visualize the h5\n",
    "trajectory with a non-hdf5 enabled app like PyMol or VMD, you can\n",
    "use mdconvert on the command line to easily transform it to NetCDF, DCD,\n",
    "or any other format of your preference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists('ala2.h5'):\n",
    "    simulation.reporters.append(mdtraj.reporters.HDF5Reporter('ala2.h5', 1000))\n",
    "    simulation.step(100000)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python",
   "language": "python",
   "name": "python"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
